A lattice study of interaction mechanisms in a heavy-light meson-meson system 
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We study mass spectra of a meson-meson system involving two light and two heavy quarks on an 
anisotropic lattice. The heavy quarks are treated in the static approximation. The dependence of the 
spectrum on the relative distance of the heavy quarks is extracted from the lattice simulation using 
the maximum entropy method (MEM). A correlation matrix of meson-meson operators emphasizing 
quark and gluon exchange degrees of freedom is employed in an attempt to learn about aspects of 
mechanisms of hadronic interaction. 
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I. INTRODUCTION 

The term "lattice hadron physics" has been coined 
for strong interaction physics based on first principles, 
i.e. quantum chromodynamics (QCD). In recent years 
hadron physics has emerged as a field in its own right 
I]. The desire to explain strong interaction phenomena 
in terms of the underlying dynamics of quarks and glu- 
ons sets the field apart from traditional nuclear physics, 
which emphasizes an effective field theory point of view. 
Among the important issues faced by hadron physics are 
baryon and meson spectroscopy, and structure, as well as 
the mechanism of their strong interaction. Those issues 
have in common the need to deal with excited states. 
In lattice QCD, which affords the most direct access to 
hadron physics at nuclear energy scales, excited states 
spectroscopy is just now moving into reach due to the 
use of anisotropic lattices, advanced analysis techniques, 
and powerful computing facilities. 

This work is concerned with learning about hadronic 
interaction mechanisms. We believe that much of the 
physics of hadronic interaction can be understood by in- 
vestigating heavy-light systems. In those the relative dis- 
tance of hadrons is a well defined quantity. The lattice 
'data' can be interpreted in terms of intuitive pictures, 
like potentials. Insight into mechanisms of the strong in- 
teraction flows from looking at excitations due to hadron- 
hadron operators, say $(t), at various relative distances 
r. Different choices for the structure of those operators in 
terms of their composition from quark and gluon fields 
may potentially point at interesting physics of the sys- 
tem, such as the importance of quark versus gluon ex- 
change degrees of freedom as a function of r. To extract 
information of this kind the computation of matrix el- 
ements (n|$(io)|0) between the vacuum |0) and ground 
and excited states \n), n > 0, will be useful. 

From the vantage point of a numerical lattice simula- 
tion this can be a notoriously difficult problem. Hadron- 
hadron operators are prone to produce very noisy corre- 
lation functions. Extracting spectral information in the 



standard fashion, i.e. trying to identify a plateau in an 
effective mass function, may not be practical. 

An alternative analysis method that is just being 'dis- 
covered' by lattice practitioners is application of the 
maximum entropy method (MEM) or an otherwise con- 
strained form of Bayesian inference |2|, ||, Qj. From the 
Bayesian perspective the parameters of a model for the 
Euclidean time correlation function are viewed as random 
variables drawn from a certain probability distribution 
function. The latter, known as the posterior probability 
V [p <— C] , is the conditional probability for a certain pa- 
rameter set p given a data set C. In our case the data 
set C is the lattice-measured time correlation function 
C(t,to) and the parameter set p is the spectral density 
function in the model 



F(p\t,t ) = / cL^He-^-'o) 



(1) 



Discretization understood, in this approach the num- 
ber of parameters is allowed to exceed the number of 
data points without causing conceptual problems. In the 
MEM the posterior probability is constructed from the 
X 2 -distance between the lattice data C and the model F, 
and the information content of p measured by the entropy 
S = — J dujp(uj) In p{u>). Usually, the result inferred from 
the Bayesian approach is the most likely parameter set p. 

Analytically, the spectral density function is a sum of 
discrete (S-peaks 



p(uj) = Y]8(u- oj n ) | (n\<f>(t 



(2) 
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From a numerical viewpoint computational constraints 
on C render the peak widths finite. Physical quantities 
are contained in each peak as low u moments. Among 
those are the peak volume |(n|$(io)|0)| 2 and the peak 
energy E n , i.e. the mean value of w. 

The Bayesian spectral analysis of the lattice data is 
an interesting problem in itself. It leads to the discus- 
sion of a number of computational strategies in a general 
context. In order to keep this presentation focused, we 
refer the reader to a separate paper || where selected 
aspects of Bayesian spectral analysis are discussed, using 
the same lattice simulation data as a testing ground. 
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Starting a little more than a decade ago lattice work 
on hadronic interaction has followed mostly two tracks; 
investigation of heavy-light systems for varying hadron 
relative distance and the computation of scattering 
phase shifts from energy spectra of two hadrons in a fi- 
nite box JjJ. We will not review these subjects here, 
but rather mention a few leads to facilitate following the 
literature.^] Extraction of scattering phase shifts has 
been successful in terms of structureless particles M . On 
the other hand, scattering of composite hadrons within 
Liischcr's framework is considerably more difficult. Scat- 
tering phase shifts have been obtained for the S-wave 
interaction of the tttt system in the 1 = 2 channel ||, 
and by the CP-PACS collaboration [[To). Considerable 
computing power was brought to bear by the JLQCD 
collaboration to extracting scattering lengths for the 7T7T, 
7rN, and NN systems O], While lattice volume limita- 
tions hinder a realistic treatment of the NN system, the it 
scattering lengths of Ref. [[n] stand out as the only quan- 
titative results for hadronic interaction from the lattice 
to date. Systems with heavy, even static, quarks will 
not yield quantitative results, but are useful to gain a 
deeper understanding of the mechanisms that lead to the 
phenomenology of nuclear forces. Green and coworkers 
have been studying energies of static four-quark systems 
in various geometric configurations, most recently with 
two static and two light quarks [|l2). Their results shed 
light on the importance of many-body versus two-body 
forces in hadrons. Other studies of heavy- light systems 
fl3f focus on adiabatic potentials, i.e. the ground state 
energies as functions of the relative hadron-hadron sep- 
aration. Along these lines, most of the hadronic inter- 
action work by members of the UKQCD collaboration 
can be traced from fl4|| . There, heavy-light meson-meson 
systems are classified according to isospin and spin con- 
figurations of the light quarks (B mesons). Depending 
on the channel, interaction potentials at around ~ 0.5fm 
come out attractive or repulsive, with magnitudes rarely 
exceeding w 50MeV. Though this is typical for nuclear 
physics, from the point of view of a lattice simulation 
these are very small energy differences to be measured. 
This is a generic problem for hadronic interaction physics 
on the lattice. 

In the present work we report on a study of a sys- 
tem of two heavy-light mesons based on two operators. 
The first one, <I>i, is the standard product of two local 
pseudoscalar heavy-light meson fields at distance r. The 
other one, $2; is similar, but nonlocal, probing color re- 
arrangement. In some sense <J>i and $2 test the impor- 
tance of quark and gluon exchange degrees of freedom for 
the interaction. They enter into a 2 x 2 time correlation 
matrix. Our goal here is to learn about the interaction 
mechanisms represented by those operators as a function 
of the relative meson-meson distance. The current work 
goes beyond previous studies of heavy- light meson- meson 
systems mainly in that (i) nonlocal operators are consid- 
ered and (ii) operator mixing is built into time correlation 
functions, allowing the system to 'dynamically pick' the 



dominant excitation mechanism as the meson-meson sep- 
aration r varies. Moreover, due to the spectral analysis 
used in this work we are able to (iii) extract the actual 
excitation strengths of the ground and excited two-meson 
states. 

Preliminary results reported in [l6) were based on 
an analysis of the diagonal correlator elements, with no 
mixing. There, averaging over annealing start configura- 
tions had not been done, the spectral density functions 
came from single annealing runs. As it turns out this is 
a source of systematic error that can not be tolerated in 
the light of the smallness of the extracted energy shifts. 



II. LATTICE ACTION 

The meson-meson operators employed in this simula- 
tion lead to somewhat massive states. The resulting steep 
drop of time correlation functions, particularly for the ex- 
cited states, makes it very difficult to analyze the lattice 
signal. To deal with this situation we use an anisotropic 
lattice action. If the aspect ratio £ = a s /at of the lattice 
constants in space and time directions, respectively, is 
made larger than one the number of usable time correla- 
tion function data is increased before the signal 'vanishes' 
into noise. Anisotropic lattices have been essential for 
computing the glue ball mass spectrum jL7j. The cur- 
rent simulation of hadronic interaction has in common 
the need for extracting excited states. 

The gauge field part of the lattice action has the form 

S G [U] = f3j2 c ^i with (3a) 

l 

{l e = ^ReTr[l - U(C)} . (3b) 

CeS e 

Here /3 = 6/g 2 in terms of the gauge coupling g, I labels 
sets Se of closed lattice contours C and U(C) denotes the 
path ordered product of gauge field link variables along 
C. We adopt a tree- level tadpole improvement scheme 
with four classes of loops: all oriented spatial elementary 
plaquettes £ = ss, temporal elementary plaquettes £ = 
si, spatial planar rectangles £ = sss, and short temporal 
planar rectangles with two spatial and one temporal link 
I = sst. Specifically 

~12^ s nsss " 12^"*" J ' (4) 
where u s and Ut are spatial and temporal link renormal- 
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ization factors |Q , and 
n ~ = E E ^ReTrl 1 - W x 

cc 1<^i<i/<3 

[/.(x + ^^^ + ^C/t^)] (5) 
= E E ^RcTr[l - C/^x) x 

x 1</j<3 

tf 4 (x + + 4) (6) 

n «- = 2 E jReTr[l-[/ M (x)C/ AI (a; + /i) x 

I7„(a: + 2/2)E/J(a; + * + £)[/£ (a; + (a:)] (7) 
n »* = E E ^ReTrll-C/^^t/^x + A) x 

x l<fj,<3 

U 4 (x + 2(i)Ul(x + 4 + pL)Ul{x + i)u\{x)] . (8) 

This action is the same as in | }L7| , the chosen parameters 
are /3 = 2.4 and £ = 3. The spatial link renormalization is 
computed self consistently from the average spatial pla- 
quette, 

u s = (l-(n ss )/3L 3 ) lf \ (9) 

while the temporal renormalization factor is set to one, 
Ut = 1. According to |l7j this results in a spatial lattice 
constant of about a s ss 0.25fm, or a^ 1 ss 800MeV. 
The fermion action is 

S F [U, lf),ii]=^2 $fA»(xi)QfAiJ,,gBv(x, y)^ g B v {y) , (10) 
x,y 

with /, <7=flavor, A, _B=color, and fi, ;/=Dirac indices. 
The fermion matrix Q is assumed flavor diagonal 

QfA^,gBu{x, y) = Sf g Q](l >Bl/ (x, V) > ( u ) 

and has a Wilson and a clover term jl9], p0| 

Q (/) = 1 - k (s \M - c sw K) . (12) 

In detail 

M(x,y)= (13) 
1 3 

+ - 74)C^(a:)« x+ a jl) + {n+l4)Ul(y)S x y+i \ . 

The critical hopping parameter in the anisotropic case is 
k c = (2£r f + 6r s )~ 1 . The Wilson parameters are chosen 
as r s — r t — 1. Because the lattice is coarse in the 
space directions and, more importantly, because we wish 
to avoid problems with ghosts (unphysical branches in 



the lattice-quark dispersion relation we use clover 
improvement [po| only in the spatial planes, 

K(x,y) = S x .y^ 0*0 -PtAxj) ■ 

s l<fj,<v<3 

(14) 

Here o>„ = §[7^,7*/] and P^(x) = | ^=1 -ffe, ^( a; ) is 
made from four transport operators along oriented 4-link 
paths in the ji—v plane, starting and terminating at x, 
for example Pi rfiU (x) = U f _ l (x)U u (x + p)U^(x + v)Ul{x), 
see @. 

The strength coefficient is c sw = 1. At k^> = 0.0679 
we have m^/m p w 0.75 on a L 3 x T = 10 3 x 30 lattice. 
Using m c sCLt = 0.4676 from [H a crude estimate for the 
quark mass puts it within 15% above the strange mass 
scale. We have used a hybrid molecular dynamics algo- 
rithm (HMC) (22) to generate Njj — 708 quenched gauge 
configurations. 

III. OPERATORS 

At this time we wish to address the physical mecha- 
nisms responsible for the features of hadronic interaction 
rather than making quantitatively precise predictions. 
Toward this end we employ, for a two-meson system, a 
set of operators meant to excite different QCD degrees 
of freedom. For two heavy-light pseudoscalar mesons, for 
example, the operator 

=^2Sf^ s _gQ A (xt)'y 5 q A (xt) Q B {yt)l5lB{yt) 

x,y 

(15) 

is built from separable products of two color singlets, 
describing two mesons at relative distance f. In ( |T5| ) 
Q and q are the heavy and light quark fields, respec- 
tively, and A, B are color indices. Color contractions are 
done between quark fields which spatially coincide at x 
and y, respectively. So far, only local operators of this 
kind have been employed for hadronic interaction studies 
H [b| . At small values of r color coupling schemes 
involving quarks in different mesons may become dynam- 
ically possible. An operator testing excitations of that 
nature is 

= Sf ,s-y Up-,AA' (xt, yt) U P ,. B , B (xt, yt) 

Q A{xt)^q B {xt) Q B ,(yt)-y 5 q A ,(yt) . (16) 

It involves link products Up(x,y) along spatial paths P 
within a fixed time slice. The operator $2(0 interpolates 
fields that are still products of two color singlets, but 
those are now formed from quarks at different locations 
x and y. 

In order to simplify the notation we will tacitly assume 
dependence of all subsequent quantities on the relative 
distance f , but suppress r in most of the expressions be- 
low. 
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We are thus lead to computing the elements 

Ctf(Mo) = (*J(i)*i(*o)>, i,j = 1,2, (17) 
of 2 x 2 time correlation matrices, one for each r, where 



$ = $ — ($), means vacuum-subtracted operators. The 
correlation matrix C(t,to) built from the operators (|l5| ) 
and ([l6|) can be worked out by way of Wick's theorem, 
symbolically 



43211234 43211432 43213214 43213412 
QqQqqQqQ = QqQqqQqQ + QqQqqQqQ + QqQqqQqQ + QqQqqQqQ . 



(18) 



where pairs nn of numbers n = 1 . . . 4 denote contractions. This leads to the following expression 
<3y(t, t ) = 2S^,(J2 Si% ff U,, AA ',BB>(t, xy)U j; DD>,cc>(to, xy) 

H B ^> Dl y(yt,yto)H AVtCX (xt,xto)[GA'^\C'\'(yt,yto)G BV:D \(B,B ) - GA>v>,D\{yt,xto)GBv,C'\>(xt,yt )}) . (19) 



In (|T^) the symmetrized Kronecker symbol 

5i + l = S+? tffl + S- ft p, (20) 



is related to 0(3,Z) symmetry, projecting heavy-quark 
distances to absolute lengths r = \x — y\. The heavy 
anti-quark propagator is employed in the static approx- 
imation, i.e. the leading term of the hopping parameter 
expansion 



1 



~ t0 ^(l + 74)A„£/ C 



CA(xt ,xt) , 

( 21 ) 

where U (xto, xt) is the path-ordered link variable product 
along a straight line from xto to xt. Finally, depending 
on whether $i or $2 is involved in the correlator matrix 
element, the contour operators 



Ul;AA',BB'(t,Xy) 
U 2 ;AA',BB'(t,Xy) 



SabSa'B' (22) 
Up-aa' (xt,yt)U P ,. BB ,(xt,yt) 
-U p ,.aa> (xt, yt)U* P:BB , (xt, yt) (23) 



are trivial, or involve path-ordered link variable prod- 
ucts along purely spatial paths P, P' from xt to yt. On 
the present lattice we consider straight on-axis paths of 
lengths r = 1,2,3,4. 

Diagrammatic likenesses of the correlator matrix ele- 
ments ([[9]) are shown in Fig. [TJ. The simplest one, the 
graph of C\\, is the standard quark exchange diagram 
that is usually considered (|, ||, [H^, [u), while the rest 
involve exp licit gluon degrees of freedom. 

In (|l^ ) light-quark propagator elements, e.g. 
G A' v' ,D\(yt, xto), between arbitrary lattice sites xt,yto 
emerge. The source time slice to may be kept fixed, how- 
ever, all-to-all spatial propagator elements are needed. 
More precisely, after working out one of the space-site 
sums in ( |T9| ) the number of propagator columns needed 
may be minimized by using translational invariance 
of the gauge field average (...). Then, the number of 



r 



propagator columns becomes equal to the number of 
relative distances f. On the other hand, employing 
all-to-all propagator elements has the advantage of 
greater flexibility with regard to varied choices for r, 
like off- axis distances (future studies), and improved 
statistics due to space-site averaging in ( J19|) . Because 
the gauge link contour operators ( |2l"| ) and (|23|) tend to 
make the correlation matrix ( |l9| ) quite noisy, we follow 
the latter strategy. 

Random-source estimation is a proven technique for 
generating all-to-all propagators ^2, 23, 24|, |2j|. Con- 
sider the following linear equation 

'^2 l '^2 l QA^B V (xx i ,yy A )X B A u M rXi) (yy A ) = 

VV4. Bi> 

Saa^ p sR^ S ^\x)6 X4X s. (24) 

On the right-hand side RS A rx ^ denotes a vector of 
length L 3 of complex random deviates. Indices super- 
scripted with s characterize the source with respect to 
color, Dirac, and time slice. We choose sources that are 
nonzero on only one fixed time slice xf = to, whereas 
for each of the 12 color-Dirac combinations A s — 1 ... 3, 
fi S = 1 ... 4, in turn, r = 1 . . . Nr different random vec- 
tors are generated. In the simulation complex Z2 dis- 
tributed random deviates were used pq] . With Nr = 8 
this results in 96 statistically independent sources per 
gauge configuration. Writing X)<r> ^ or ^ ne ensemble av- 
erage, of which j^- y~]^—i is a truncation, we have 



R (A " 

<r> 



f)(2)#V«f>^ = f (25) 



for all A s , fi s , xf with appropriately chosen normaliza- 



tion. Application of @5p to ( |24[ ) yields 



G B v,An(yt,xt ) 



J2 x^rto) (^ )jR (^rto)* ( f ) ( 26 ) 
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y t x t 



c, = 



11 



c,„ = 



12 




with m € N, and the matrix 



c 21 = 



c 





FIG. 1: Diagrammatic representation of the elements of the 
2x2 correlator matrix ([L9|) . Thick and thin solid lines mean 
heavy and light quark propagators, respectively. Gluon prop- 
agation, i.e. products of link variables, are shown as dashed 
lines. Pairs of dots refer to the same space-time point. 



as a stochastic estimator for spatial all-to-all propagator 
matrix elements. 



We have used the BiCGStab algorithm |27], [28) for solv- 
ing (H). 

Operator smearing [p9[ is a technique for enhancing the 
amplitude of the ground state in a correlation function. 
We smear the light quark fields only, defining iteratively 



ICAB{x,y) = 8AB5x,y + Ot'^2 l \U l i,AB{xt)53$-f 1 

+ul AB (yt)^ +i x] • (28) 



The real number a and the maximum value M for iter- 
ations m — . . . M are parameters. We have also used 
APE type |Fj fuzzy link variables U £ SU(3) in @ 
with the same parameter set. The above iterative pre- 
scription translates directly to the random source and 
solution vectors R and X, i.e. replacing R — > R^ M ^ 
and X X* M > in (||) yields the propagator G* M > for 
smeared fermion fields. For more technical details see pi. 

Smearing was employed as a matter of course. In fact 
the operator ([l6]) is designed with excited states of the 
two- meson system in mind. Thus light-quark field smear- 
ing and link variable fuzzing may be of limited value for 
the task at hand. We have thus chosen somewhat con- 
servative parameter values a = 1.4 and M — 4 for the 
current simulation. 



IV. SPECTRAL DENSITY ANALYSIS 

We now turn to the problem of extracting spectral in- 
formation from the correlation matrix ( ]l7| , ^9| ) . It has the 
decomposition 



"3 



cj> ni = (0|$l(t o )|n), 



B n (t-t ) 



< - (29) 



(30) 



where \n) denotes a complete set of states with energies 
E ni some of them may be negative in a lattice simulation. 
The sum in ( |2g| ) is truncated in practice, n = 1 . . . N. 
The truncation N is determined by the physics of the 
system (ultimately the lattice action) , the simulation pa- 
rameters (most importantly the lattice energy cutoff) , 
and by the choice of operators. The matrix elements (pfl) 
can be interpreted as components of N Af -dimensional 
vectors 4> n , where M is the size of the correlation matrix 
C(t,to)- Following Liischer and Wolff we diagonalize 
the correlation matrix separately on each time slice, say 

M 

}^ Cjj (t, t )v m j(t, t ) = C m (t,to)v mi (t,t ) (31a) 
i=i 

A I 

'^2 v mi{ t i t a) v m'i{t,t<d) = 5 mm > , (31b) 



q A °\xt) = q A (x) (27a) 
q A m \xt) = J2T,' C ^y)q B m - 1} m, (27b) 



denoting the eigenvectors by v m and the eigenvalues by 
C m , m = 1 . . . M. Now assume the following: 

1. The energies E n in ( p9| ) are non degenerate, and 
ordered E\ <E 2 < ... < E N . 
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2. The vectors <f> n are linearly independent, this im- 
plies N < M. 

3. There is a tc such that for all t>tc the eigenvalues 
are ordered C\ > C% > ■ ■ ■ > Cm- 

Under those conditions a theorem proven in |3lj states 
that for all n = 1 ... TV 

lim C n (t,t ) = Z^-^-^ll + Oie-^^-^)] , (32) 

t — >OG 

where Z n > and AE n = min„/j n {|.E n / — E n \} is the 
distance to the energy closest to E n . We are interested in 
the structure of the spectral rep resen tatio n of C n {t,to). 
Toward this end, applying ( plb| ) to ( 31a ) and then in- 
serting (p9|), one obtains 

M 

C m (t, t ) = J2\J2 M<M 2 e- £ " (t -' o) ■ (33) 

n^O i=l 

In addition to items 1.-3. above we will also assume: 

4. In the large -t limit the first N' of the eigenvectors 
of (31a, 31b ) converge in the sense that in 



lim v mi (t,t ) = ri m (t,t ) v ml , m= 1...N' , (34) 

t — >oo 

the vectors v m are constant, being multiplied by a 
i-dependent phase factor, \r] m (t,to)\ = 1. 



Note that N' < N < M. Thus, using (g) and @) in 
(^) we arrive at 

M 

E lE<^~| 2e " B " ( " t0) = Z m e' E -^ , (35) 

for m = 1 ... N'. Since all E n are different the exponen- 
tials are linearly independent functions of t, hence 



A I 



i E w ™^ ni i 2 = ZnSr > 



m = 1 . . . N' . 



(36) 



i=i 

The square root of this is 

M 



E< 



CrnV Z m 8 mn , m — 1 . . . N . 



(37) 



i=l 



with \Q m \ = 1. The eigenvectors of fl31a ,31b) satisfy a 
completeness relation in M-dimensional space. We split 
it into two parts 



N' 



A I 



E V mj(t,t )v; ni (t,t ) + E U fci(*!*oKi(*,*o) = Sji 



m—1 



k>N' 



(38) 

According to the assumption (|3^) all terms in the first 
sum will individually converge in the large-i limit. The 



individual terms in the second sum might not, however, 
it must of course become ^-independent as a whole, 



A I 



%= lim E v kA t ^o)v* kl {t,t„) . 



(39) 



k>N' 



Clearly the projector II = II 2 = IT is orthogonal on the 
space defined by the span of the v m , 



A I 



En. 



= 0, for m = 1 . . . N' . 



(40) 



Thus, as t — > oo, equation ( pq) assumes the form 



(41) 



N' 



on both sides of 



AT' 

^ ] v mjV ri 

771=1 

Finally, operating with 
( p7| ) and then using ( f4l| ) gives 

K = ( n VZ~v n +H(f, n , for n = l...N'. (42) 

This result relates the matrix elements <p n i — (0|$^(io)|n) 
to the solutions of the i-dependent eigenvalue problem 
ila|. 31b) in the large-t limit. An immediate conse- 



quence, derived by taking the square of ([4 2D, is (fy n (j)n 
Z n + <j>lll<frn, or 



A I 



El^l*i(*o)|0)| 2 = Z„ + ||n^„|| 2 , for n=l...N'. 

i=l 

(43) 

Thus Z n is a lower bound on the total probability for (in- 
coherent) excitations by a set of operators $j, i = 1 . . . M, 
into a certain state \n). In principle the value of N' 
can be computed from (|i|). In practice this is hard to 
accomplish, because the components of the eigenvectors 
fluctuate strongly. If the set of M operators couples to 
all N available physical states it seems reasonable to ex- 
pect that N' — N . In case that N < M we have used 
more operators than quantum states are available in the 
system. Then the projector term compensates for over 
counting the physical degrees of freedom. 

Aside from ( ^3[ ) there is an alternative way of inter- 
preting the Z n . Based on |34|) define the meson- meson 
fields 



A I 



*m(*) = E v mi®i{t) for 771 = 1 
i=l 

and consider the N' x N' correlator matrix 

Anm'(Mo) = (*LW*m'(*0)>- 



,N'. 



Inserting (44) into (|45|) , then using (|29j) and (|E 
straightforward to show that 



Dn 



'(Mo) 



,Z p--Bm(t-to) 



J mra' ^ra 



(44) 

(45) 
it is 

(46) 



On the other hand, starting from (fig), the diagonal ele- 
ments have the standard decomposition 

D mm (t,to) = E | <«| * TO (*o) |0> |a e -«- C*-*o) . (47) 

Comparison of ( f46| ) and (f47|), using linear independence 
of the exponentials again, then gives 

\(n\V rn (t Q )\0)\ 2 =5 nm Z n for n,m = 1 . . (48) 



Thus Z n = \(n\V n (h 



is the excitation probability 



of the state \n) due to an operator ^ n (t) that is optimal, 
with regard to \n), within the linear space of the original 
set $i(t),i= 1...M. 



V. DISCUSSION OF RESULTS 

The significance of the above is that it suggests an anal- 
ysis strategy for the spectral features of the two-meson 
system: Diagonalize the time correlation matrix ( |l9| ) on 
each time slice, 



Cii(t)to) Ci2(t, to) 
C2l(t>*o) C22\t, to) 



C m (t, to) 



v m i(t, to 

v m2{t, to 

Uml(Mo) 
Wm2(*,*o) 



(49) 



m = 1,2. Then view each eigenvalue as a separate time 
correlation function C m (t, to) subject to spectral analy- 
sis. In particular, seek to extract a spectral representa- 
tion of the form J du; p(w) e -w(*-*o) see ^ This can 
be done by Bayesian inference ]s|, [33| . The expected 
structure of the spectral density p is a linear combina- 
tion of peaks, see (^). In an actual numerical simulation 
those will have finite widths. More importantly, because 
the main results of Sec. IV, including ( |43| ) and (f48|), re- 
quire the limit t — > oo only the lowest-energy peak of 
each C rn (t,to) is significant. Physical information that 
can be extracted from the peak includes the energy E n 
of the state \n), and the strength Z n of excitations by 
means of the set of operators employed. 

The eigenvalue correlators C m (t,to),m = 1,2, of ( fl9| ) 
are displayed in Fig. ^ for four meson-meson relative 
distances r = 1 ... 4. Due to the non-local structure 
of the two-meson operators, leading to loop-loop cor- 
relations, somewhat noisy data are unavoidable. Link 
variable fuzzing and operator smearing was used to be 
able to use 'earlier' time slices, recall Sect. [II. Not sur- 



prisingly, the lower eigenvalue C2(t,to) is more suscepti- 
ble to noise than C\(t, to). Backward going propagation 
is present, though suppressed by four to five orders of 
magnitude. j43j By inspection of Fig. || it is apparent that 
the condition Ci(t,to) > C*2(£, £o) is fulfilled for t > to- 
This is directly obvious for most of the forward temporal 
range, say < t < 20, and for 20 < t < 30 it can be in- 
ferred from the global fits to Ci(t, to). Indicative of a vi- 
olation would be the possibility of being able to smoothly 




FIG. 2: Eigenvalue correlation functions at meson-meson rel- 
ative distances r = 1, 2, 3, 4. We show C\ (t, to) as filled circles, 
and Ci(t, to) as filled squares. Indicated are statistical errors 
from Nu — 708 gauge configurations. A missing plot symbol 
means that the statistical error exceeds the value of a data 
point. The lines are plots of the model ( |53| using the Bayesian 
results for the spectral densities p from Figs. |^and ^ 



connect two sets of consecutive data points in such a way 
that the two resulting curves would intersect. |Q We thus 
observe that the above ordering of eigenvalues, as stated 
in Sec. IV, is satisfied. 



For asymptotic times individual components of the 
eigenvectors in ( fl9|) show large statistical fluctuations. 
In Fig. U a typical example of the (complex) values of 
z = v n i(t,to)v^ 2 (t,to) is displayed. This quantity is rel- 
evant to the assumption (|34|). In the region of forward 
propagation, say t < 20, the real part of z exhibits in- 
creasing stochastic errors as t becomes large, but sta- 
bilizes within those bounds. A similar statement can be 
made for the imaginary part of z, adding that it is consis- 
tent with zero. This justifies ( p4[) w ithin the limitations 
given by the quality of the data. Ea| 

In Fig. H we show the magnitudes |«n(i, £o)| 2 of eigen- 
vector components versus t for four relative distances. 
Substantial noise notwithstanding we have strictly ap- 
plied the prescription of Sect. IV and computed the eigen- 
value correlators following ( 314311 ) verbatim 



C m (t, to) 



^ Vmi(t,to)Cij(t,t )v m j(t,to) , (50) 



8 




10 15 20 25 30 

t. 



FIG. 3: Real and imaginary parts of z = v* 1 (t,to)v„2(t,to), 
the example is for the ground state n = 1 and relative distance 
r = 2. The numerical constancy of z within errors (excluding 
the backward propagation in the region t > 20) illustrates the 
validity of (g). 



for to = 1,2. The corresponding spectral analysis is dis- 
cussed below. 

Stochastic fluctuations of the eigenvector components 
are enhanced by diagonalizing (|49|) separately on each 
time slice. An attractive alternative is to replace the 
eigenvector components v m i(t, to) in ( |so|) with time aver- 
aged components v m i taken at asymptotic times. Within 
the present data set the available time window is 10 < 
t < 20, excluding backward propagation. We have made 
fits to |fii (t, to) 1 2 with a constant model s in the time 
slice range 11 < t < 17, and also with a control set in 
11 < t < 19. The results are listed in Tab. | Specific 
values for s give rise to, time independent, vectors 



s/T 



V 2 



(51) 



up to an arbitrary phase. Thus we have also performed 
a spectral analysis based on the correlators 



C m (t, to) 



2 

£ 



(52) 



for m = 1,2. Besides smoothing out statistical fluctu- 
ations of the eigenvector components the advantage of 
(|52|) is that the asymptotic form of the correlator J5C| ) 
is now used on all time slices. This should improve the 
signal derived from the Bayesian spectral analysis which 
makes use of data on all times slices. 

The values of s inform us about operator mixing as 
the relative distance r changes. They are a measure of 
how strongly the operator $i couples to a meson-meson 
system in the ground state. As Fig. |s| shows this mea- 
sure distinctly decreases from about 1.0 as r becomes 
smaller. Since $i is designed to test quark exchange de- 
grees of freedom we see that those gradually become less 
important at smaller relative distances. By the same to- 
ken 1 — s measures the ground state coupling strength of 
$2- Thus we learn that quark exchange is the dominant 
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FIG. 4: Time dependence of the squared magnitudes y = 
\v m i(t, to ) | 2 of the 7 = 1 component of the ground state m = 1, 
for relative distances r = 1, 2, 3, 4. The uncertainties are jack- 
knife standard errors. The horizontal lines indicate constant 
model fits in the range 11 < t < 17, s from Tab [E| 



TABLE I: Time averages s of |t>n(i, to)\ 2 and the correspond- 
ing variances As. Fits for two time slice ranges are listed. 





11 < t < 17 


11 < t < 19 


r 


s 


As 


s As 


1.0 


0.629 


0.132 


0.645 0.138 


2.0 


0.854 


0.060 


0.831 0.134 


3.0 


0.913 


0.173 


0.861 0.259 


4.0 


0.996 


0.033 


0.994 0.046 



interaction mechanism at large r, while gluon exchange 
gradually takes over as r decreases. A glance at Fig. B re- 
veals that the mechanisms become balanced (s w 0.5) at 
distances r somewhat less than 0.5, or 0.25fm in physical 
units. 

The Bayesian analysis of time correlation functions has 
been extensively discussed in g, using the same lattice 
(raw) data. We here briefly state the main points for 
coherence of presentation, but otherwise refer the reader 
to H . We expect the lattice data to fit the model (Q) , or 
rather its discretized form 



F(p\t,t ) 



E P ke 

k=K- 



-u k (t—t ) 



(53) 
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FIG. 5: Asymptotic time averaged values s of \vn(t,to)\ 2 , 
and their variances, as functions of the relative meson-meson 
distance r. Results for two time slice ranges are shown. The 
dashed curves are quadratic polynomial fits. 



with pk = Aojp(u>k) and ojk = Acjfc. pq] The objective is 
to compute the spectral density function p. Toward this 
end consider the functional 



W[p] = \ X 2 [p\ 



aS[p] . 



(54) 



The first term involves the usual x 2 -distance between the 
lattice data C m (t, to) and the model F(p\t, to), computed 
with the full covariance matrix derived from gauge field 
configuration statistics, and 



K 4 



k=K- 



S[p] = ^2 [Pk - m k - Pk In 



Pk 



(55) 



is the (Shannon- Jaynes) entropy |3J, [55]]. It plays the role 
of a Bayesian prior [32]. The configuration m = {mk ■ 
K- < k < K±} is called the default model. Another pa- 
rameter in (p4) is the entropy strength a. The optimiza- 
tion problem x 2 [p] — min is ambiguous because in prac- 
tice a reasonable resolution Aoj will result in the number 
K _|_ — K_ + 1 of fit parameters pk being much larger than 
the number T of simulation data. However, W[p] has a 
unique absolute minimum From the viewpoint of 

Bayesian statistics p is interpreted as a random variable 
subject to a certain probability distribution (posterior 
probability) . The most likely p is the one that minimizes 
W^[p]. Finding the minimum within this framework is 
known as the maximum entropy method (MEM). It is 
designed to minimize the information not supported by 
the data. Loosely speaking, we seek to minimize X 2 [p] 
while assuming as little information as possible about 
the spectral density p. To solve the optimization prob- 
lem W[p] — min, probabilistic methods seem closest in 
spirit to the Bayesian stochastic interpretation of p. We 
have thus employed simulated annealing p6[ , or cooling, 
based on the partition function 



Z w = J [dp]e 



-0wW[p] 



In g] we have studied the dependence of the resulting 
spectral density on (i) the entropy strength parameter 
a, (ii) the default model m, and (iii) the annealing start 
configuration. It was found that p was essentially inde- 
pendent of a and of to, and that the expectation values of 
low oj moments were independent of the annealing start 
within errors native to the lattice data set. 

Results of the MEM analysis are presented in Fig. ^ 
for the eigenvalue correlators ([30]) and in Fig. ^| for the 
asymptotic stabilized correlators (|5^). As discussed in 
Q the parameter choices are not critical. Specifically, 
the entropy strength is a = 5.0 x 10~ 7 , and the default 
model is constant with mk = 10 -12 , < k < K + . The 
annealing schedule is given in Jj|. All spectral densities 
are averages over eight random annealing starts. The 
uj discretization is set by Aui = 0.02, and K- = —100, 
K + = +200. Those numbers reflect the lattice design, 
like the energy cutoff at -1 , and other considerations, see 
Q. The oj interval is larger and the resolution much 
finer compared to the discretization used in jjjj. Note 
that the entire spectral mass range —2.0 < u> < +4.0, 
including backward going propagation, is utilized in the 
spectral analysis. This is true of both the eigenvalue cor- 
relators and the asymptotic stabilized correlators. Most 
of the spectral structure,however, is invisible on the lin- 
ear scales used in Figs, jfi] and 0, particularly for u> < 
(backward propagation) . 

As discussed in Sect. |y], in the limit t — > oo only 
the lowest oj peak from each correlator C m (t,to),m = 
1 . . . M , should be used to extract physical information. 
We will refer to those as primary peaks. Suppose pri- 
mary peaks are seen in the spectral densities belonging 
to C n (t, to), for n = 1 . . . N < M. We loosely character- 
ize those by S n — {oj : oj £ peak #n}. In ^ it was argued 
that low oj moments of the spectral densities p(oj) can be 
reliably extracted. Specifically, these are the peak vol- 
ume Z n , the mean energy E n , and the width A n of the 
peak 



A 2 



dojp(oj) 



Z n 1 J s d0jp{0j)0J 



= z- 



<5n 



dojpioj) {oj — E n ) 



(57) 
(58) 
(59) 



(56) 



Note that the theorem (32) establishes the peak vol- 
ume rt57 ] ) as identical with the factors Z n introduced 
in Sect. [TV|, the caveat being that numerically extracted 
peaks have finite widths. In particular, the peak volumes 
(|S7]) have the interpretation given by (fl3]), or by ([i"8| ) as 
optimal excitation probabilities. 

In Figs. |^ and 7] the (lowest mass) primary peaks 
clearly dominate both ground state spectral density func- 
tions, C\ and C\, for each r = 1 . . .4. With reference to 
( |43| ) the secondary peaks at larger mass may indicate 
that ||II0x|| 2 > 0, their volumes are smaller though. It is 
much harder to make out a distinct peak structure in the 
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FIG. 6: Spectral density functions p of the eigenvalue correla- 
tor dic| ) obtained by way of simulated annealing. The graphs 
are the average over eight random annealing starts. Spec- 
tra are shown side-by-side for the ground state correlator Ci 
and the excited state correlator C2 for meson-meson relative 
distances r = 1, 2, 3, 4. 



FIG. 7: Spectral density functions as in Fig, Q, but for the 
asymptotic stabilized correlator functions (p2j). 
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spectral density functions for the excited state eigenvalue 
correlators C 2 . We attribute this fact to strong statisti- 
cal fluctuations of the eigenvector components v m i(t, to), 
spoiling the signal. The asymptotic stabilized excited 
state correlators Ci alleviate that problem, see Fig. |?|. 
The spectral peaks of C2 are broad, for r > 1. This is an 
indication that the lattice data support them with only 
a few consecutive points in the time correlation func- 
tion. In other words, there is not enough information 
in the data for distinct narrow peaks to develop against 
the entropy background. Although the peaks are wider, 
lower mass primary peaks are clearly distinguishable for 
r = 1 ... 4 in Fig. [| In Tabs. |l| and |FFi] we list the vol- 
umes, energies, and widths of all primary peaks of Figs. ^ 
and 0- For the excited states the C2 data do not always 
clearly define a primary peak. The corresponding num- 
bers in Tab. use w cuts of 1.32, 1.58, 3.20 and 2.80, for 
r = 1,2,3 and 4, respectively. On the other hand, the C2 
data originating with the asymptotic stabilized correla- 
tors provide us with a much improved picture. There, by 
inspection of Fig. 0, the ui cuts are 1.28, 1.72, 2.20 and 
2.46, for r = 1 ... 4, respectively. Most of the subsequent 
discussion therefore uses results of the C analysis. 

With an appropriate annealing schedule cooling fluctu- 
ations at the final 'temperature' can be made negligible. 
The size of the statistical error from Nu — 708 gauge 
configurations is comparable to uncertainties from the 
MEM analysis. This was looked at in ||. Those un- 
certainties are introduced through using different start 
configurations p in the annealing process. There is no 
magic way to eliminate those, nor would it be desirable, 
because they test a property of the data set. Although 
theoretically W[p] has a unique absolute minimum sole 
knowledge of its location, say p m i n , is deficient. It should 
be supplemented by having some notion about the shal- 
lowness, or as the case may be, the distribution of local 
minima close to p m i n with values W[p] not much different 
from VFfpmin] which the annealing algorithm may settle 
into. Those manifest themselves in micro fluctuations 
(on a scale of Aw) of the spectral density functions. Lo- 
cal minima close to p ml!l appear to be numerous. The 



numbers in Tabs. |Tj and III are from averages over eight 
annealing start configurations, and the uncertainties are 
the corresponding standard deviations. They are indica- 
tive of the spread of local minima of W[p] in the vicinity 
of Pmin and, ultimately of the uncertainty of the results. 

In view of the broader peaks in Figs. ^ and [7] it may be 
argued that the widths A2 and A2 supersede the stan- 
dard deviations derived from annealing start configura- 
tions as a useful indicator for the uncertainties of the 
average energies E2 and E2, see Tabs. || and H. In 
fact there are at least three types of indicators: (i) the 
gauge field configuration statistical error, (ii) the anneal- 
ing start configuration standard deviation, and (iii) the 
peak width. They all point at different aspects of the 
uncertainty of the E n . For example, the A n convey the 
aspect of information content of the data in the sense of 



TABLE II: Low ui moment features of the primary 
spectral peaks extracted from the eigenvalue correlators 
C m (t,to),m = 1,2 for different relative meson-meson dis- 
tances r. The corresponding spectral density functions are 
displayed in Fig. ^| Listed are the peak volume Z m , the peak 
energy Em, and the peak width A m , for m = 1,2, as de- 



fined in (|57j)-(|5E|). All entries are averages over eight random 
annealing start configurations, the uncertainties are the cor- 
responding standard deviations. 



Zi 



Ai 



Z 2 



E 2 



A, 



1.0 124.3(5) 0.985(2) 0.062(5) 97.5(1) 1.277(2) 0.263(3) 

2.0 141.(2) 0.947(5) 0.066(9) 129.0(4) 1.851(6) 0.595(9) 

3.0 176.(2) 0.971(4) 0.083(9) 119.7(5) 2.109(8) 0.569(7) 

4.0 174.(1) 0.970(4) 0.093(6) 84.7(1) 1.983(2) 0.203(3) 



TABLE III: Low uj moment features of the primary spectral 
peaks like Tab. [| but for the asymptotic stabilized correla- 
tors C m (t,to),m = 1, 2 of ([52]). The corresponding spectral 
density functions are displayed in Fig. [?]. 



r 


Zx 


Ei 


Ax 


z 2 


E 2 


A 2 


1.0 


107.(1) 


0.997(3) 


0.072(5) 


67.2(9) 


1.002(5) 


0.072(9) 


2.0 


141.(2) 


0.951(4) 


0.074(7) 


44.(1) 


1.29(1) 


0.191(8) 


3.0 


175.7(9) 


0.974(3) 


0.090(7) 


60.(2) 


1.63(1) 


0.246(9) 


4.0 


162.6(6) 


0.969(2) 


0.094(3) 


65.(1) 


1.64(1) 


0.263(9) 



[B4| . A deeper discussion is beyond the scope of this work. 
For the current results the gauge field configuration sta- 
tistical error is typically smaller than 2% and thus much 
less than the other two measures of uncertainty. Then, 
the annealing start configuration standard deviation in 
most cases is much less than the peak width. We here 
adopt the point of view that, since Bayesian inference is 
built on the information content of the data, the peak 
width is the appropriate measure of uncertainty. How- 
ever, for reference to (i) refer to S, and we include (ii) 
in the results below as appropriate. 

The ground state energies E\ of the meson-meson sys- 
tem are plotted in Fig. g| together with the mass 2m and 
the error band for non-interacting mesons. The mass 
to = 0.4676 and the gauge configuration statistical error 
Am = 0.0075 are those from m e ff,o injj- Error bars on 
the data points are the peak widths Ai. Although the 
errors extend well below V — 0, a systematic tendency 
for repulsion, averaging about w 90MeV above 2m, is 
apparent. These features are consistent with previous 
lattice studies of the pseudoscalar meson-meson system 
in the same isospin channel [^| . 

We also display in Fig. || results from |Q which relate 
to the current discussion. In Jl4| operators are charac- 
terized by the isospin and spin symmetry of the light 
quarks. Because in the present work all the light quarks 
have the same flavor, only the symmetric isospin combi- 
nation is relevant; in the notation of Jl4| this is I = 1. 
Figure || shows their S = and S = 1 results from two 
lattices. The uncertainties are gauge configuration sta- 
tistical errors and are thus much smaller than our peak 
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FIG. 8: Ground state energies of the meson-meson system 
for various relative distances r. The mass 2m of the non- 
interacting system and the gauge configuration statistical er- 
ror, both from Q|, are shown as horizontal lines. The solid 
plot symbols correspond to the energies E 1 of this work, the 
uncertainties are spectral peak widths Ai. For comparison 
we show the 7 = 1 and S — 0, 1 data for two lattices compiled 
from 0. 



width error bars. The energy values populate a ±50MeV 
band around V = 0. It should be noted that the au- 
thors of jl4j go through extraordinary length to obtain 
the 'best' ground state energies possible by means of di- 
agonalizing a matrix of operators (similar to $i) with 
different fuzzing levels of the gauge fields. This goes be- 
yond the widely used practice to employ a few iterative 
steps, say N, of gauge link fuzzing Q and quark field 
smearing p{| . A reasonable choice for the number of it- 
erations is to make Na s about the radius of the hadron 
considered. (In our case, with a s — 0.25fm, this would 
be N — 4, for example.) The result is a spatially ex- 
tended operator, spreading across a hadron volume, that 
makes correlation functions assume asymptotic behavior 
at 'earlier' time slices. This is the procedure adopted in 
the current work. This practice also has the side effect of 
lowering the ground state energy obtained from the simu- 
lation because, numerically, contaminations from excited 
states are reduced. This effect is of course enhanced if 
the ground state is extracted from diagonalizing a corre- 
lation matrix from operators of different fuzzing levels, 
as employed in fllij . The size of the effect can be seen 
in Fig. ||. We note that the physics goal of the authors 
of nM is to investigate if the BB system can possibly be 
bound. In that context even small effects on the ground 
state energy are vital, so the elaborate matrix fuzzing 
procedure is justified. However, the physics goals of this 
work are completely different. Studying the interaction 
mechanism rests on a comparison of ground and excited 
state energies, as r changes. Because the excited states 
energy levels are substantially larger, tiny shifts in the 
ground state energies are irrelevant. 

This situation is reflected in Fig. where we show 



FIG. 9: Energies of the excited (squares) and ground (circles) 
meson-meson states relative to twice the single meson mass 
in physical units. The uncertainties are peak widths from the 
spectral analysis. The data points at r — 0.25fm are shifted 
sideways slightly, for better visibility. The curve is a fit with 
y = —a/r + br + c, see text. 



excited state and ground state energies together, ver- 
sus the relative meson-meson distance. The energy E2 
drops considerably as r decreases. At large r, accord- 
ing to Tab. U and Fig. ||, the operator 5>2(t) is mainly 
responsible for the values of Ei- Bearing in mind that 
<I>2(t), see ( M ) and Fig. |l involves color charges sep- 
arated by a distance r, which should be confined, we 
have fitted the data with the model y = —a/r + br + c, 
where the parameter b was fixed to the string tension, 
Vb = 0.44GeV or b = 0.968GeV/fm. The remaining fit 
returns a = 0.30(14) GeVfm, c = 1.09(52) GeV, where 
the uncertainties are variances. The resulting curve is 
shown in Fig. ^|. Around r rs 0.2fm a level crossing be- 
tween the ground and excited state energies of the meson- 
meson system apparently occurs. From the viewpoint of 
an adiabatic potential this distance defines the transi- 
tion point (avoided level crossing) between a weakly re- 
pulsive interaction mostly mediated by quark exchange, 
and strongly attractive interaction from gluon exchange 
degrees of freedom. The transition distance is consistent 
with Fig. |^, where at a point somewhat smaller than 
r ~ 0.5 marks equal s values for the $1 and $2 com- 
ponents of the asymptotic eigenvectors v\ and Vi- The 
picture emerging from the current results thus is that of 
heavy-light pseudoscalar mesons having a radius of about 
R w O.lfm, with a sharp boundary, as probed by their 
mutual interaction. Note that the dramatic change of the 
adiabatic potential, at the distance of the avoided level 
crossing nevertheless comes from a smooth transition be- 
tween the interaction mechanisms, as manifest in Fig. |^, 

By extrapolation it appears that below r ~ 0.2fm the 
interaction turns attractive. This observation is in line 
with |l4| where strong attraction at r = due to gluonic 
effects is observed, and also with results in |)| [ll|. A sim- 
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pie explanation is suggested by the SU(3) color content 
of the two-body operators. Using standard nomenclature 
p7|, B8[ we note that only the singlet from 



3 = 8< 



(60) 



is used in the construction of JlE|). Thus, generically, the 
color-source structure of the operator <I>i is 



(61) 



where 1 and 2 denote the two color sources. The decom- 
position of the product of two color octets 



8<g>8 = 27 10 08080 10 01 



(62) 



also contains a singlet and therefore mixes with ( J6l| ) . The 
construction of (hq) involves two gauge link paths. Thus 
in some sense the color-source structure of $2 is, schemat- 
ically, described by the Clebsch-Gordon series 



$^(1,2) 



E 

hi 



8 8 

i j 



4 8) (l)0f (2). (63) 



The interaction energies for one-gluon exchange in those 
states are proportional to the expectation values of 
2F(1) • F{2) = (F(l) + F(2)) 2 - F(l) 2 ~ F(2) 2 where 
F 2 is an SU(3) Casimir operator |Q. A simple calcula- 
tion gives 



(2F(1)-F(2))^ = 
(2F(1)-F(2)) $(1) = -6. 



(64) 
(65) 



which indicates (possibly strong) attractive interaction 
at small r such as described by excitations due to $2. 
This is only a schematic picture, but is it appears to be 
consistent with the results of the lattice simulation. [ft7| 

Finally, we turn to the transition matrix elements. 
Equations (^) and ( p8| ) relate the peak volumes Z n to 
(unnormalized) excitation probabilities of the states \n), 



Z n = |<nKl$l(t )+«„2$2(*o)|0)| ; 

< |<»|* 1 (to)|0>| 2 + |<n|* a (*o)|0>| 



(66a) 
(66b) 



The coefficients v n i in ( p6a[ ) , being eigenvector compo- 
nents, ensure that Z n is maximal within the linear space 
of the available operators $i(to),i — 1,2, see Sect. 0. 
Accordingly we interpret Z n , or 



Zn = Z n /{Zi + Z 2 ) ■ 



(67) 



as measures for the effectiveness of the set of operators to 
excite the state \n). The corresponding results, compiled 
from Tab. Ill, are displayed in Fig. [It]. The decrease of 
Z\ as r becomes smaller means that the operators $1 and 
$2 start failing to capture some physics of the two-meson 
ground state at small relative distance. This is a sign of 
an emerging interaction mechanism that is is not well 




1 1 0.0 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 
r [to] r [fm] 

FIG. 10: Dependence of the peak volumes Z n and the nor- 
malized peak volumes z n for the ground and excited states 
n = 1,2, The uncertainties, which are annealing start stan- 
dard deviations, are small and obscured by the plotting sym- 
bols. 



represented by $1,2- Other operators should eventually 
be added to the set. From Fig. O we also see that the 
operator set $1,2 is about 60-70% effective in creating 
the two-meson ground state, and the rest is left to ex- 
cited state creation. It is possible that gauge link fuzzing 
and quark field smearing is responsible for the enhanced 
ground state presence. The normalized excitation rates 
are essentially independent of r, as Fig. |10| shows. 



VI. SUMMARY AND CONCLUSION 

In an attempt to learn about the physics of hadronic in- 
teraction we have studied the spectrum of a meson-meson 
system consisting of two light and two static quarks. The 
static quarks define the relative distance r between the 
mesons. We have supplemented the commonly used lo- 
cal two-meson operator $i(t) with a nonlocal one <t>2 (i) 
which is built from spatially extended color singlets. The 
eigenvalues of the corresponding 2x2 time correlation 
matrix were used in the spectral analysis thus allowing 
operator mixing. We have employed the maximum en- 
tropy method (MEM), a form of Bayesian inference, to 
extract the spectral densities of the ground and the ex- 
cited states from the time correlation matrix. 

The spectral analysis yields both energies of the ground 
and excited states of the meson-meson system as a func- 
tion of the relative distance r. While, at large r, the 
ground state energy is weakly repulsive and flat, the ex- 
cited state level is strongly r dependent and decreases 
substantially as r becomes small. By way of extrapola- 
tion, the salient feature of the spectrum is a level cross- 
ing of the ground and excited states at about r w 0.2fm. 
There, the adiabatic potential changes from weak repul- 
sion to strong attraction. Analysis of the eigenvalues of 
the correlation matrix, at asymptotic times, reveals that 
the interaction mechanism gradually changes from being 
dominated by quark degrees of freedom (quark-antiquark 
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exchange) at large r to gluon degrees of freedom( gluon 
pair exchange) as r becomes smaller. This view is based 
on monitoring $1.2 operator mixing as a function of r. 

A side aspect of the current work relates to the var- 
ious types of errors emerging in the analysis procedure. 
Between the gauge configuration statistical errors (which 
are small), the uncertainties native to the MEM analysis 
(which are comparable), and the widths of the spectral 
peaks (which are typically large) it remains a matter of 
judgment to decide which type of uncertainty is physi- 
cally relevant. We have here taken the point of view that 
the spectral peak width should be considered the error 
of the simulation results because it is a measure of the 
information content in the spirit the Shannon- Jaynes en- 
tropy contained in the lattice data. By this measure we 
accept that the uncertainties generally exceed the gauge 
configuration statistical errors. 

From the peak volumes of the spectral density func- 
tions we learn about the efficacy of the operators $12 
of coupling to the low-lying physical excitations of the 



meson-meson system. Although most of the important 
excitation mechanisms appear to be covered by $ 12 more 
operators should be employed in more detailed studies of 
hadronic interaction mechanisms. 

Technically, investigations into hadronic interaction 
are very demanding for both the need to extract small 
energy differences and having to dealing with (noisy) 
nonlocal operators involving loops that are a combina- 
tion of gauge link products and light quark propaga- 
tors. Anisotropic lattices and advanced analysis tech- 
niques seem essential tools in future studies. 
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